Study 1: Extra (not reported)

Luhrmann, Weisman, et al.

In this notebook, we conduct additional analyses for Study 1: In-depth interviews with people of faith (not reported in the main text or the Supplemental Information of the paper).

source("../../scripts_general/dependencies.R")
source("../../scripts_general/custom_funs.R")
source("../../scripts_general/var_recode_contrast.R")
source("../scripts_s1/s1_var_groups.R")
source("../../scripts_general/data_load.R")

Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  study = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  quad = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_school = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  godexpviaawe_freq = col_character(),
  sleephabit = col_character(),
  researcher_date = col_date(format = ""),
  researcher_notes = col_character(),
  site = col_character(),
  religion = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ctry = col_character(),
  wher = col_character(),
  recr = col_character(),
  whoc = col_character(),
  demo_chur = col_character(),
  demo_ethn = col_character(),
  demo_maj = col_character(),
  demo_pocc = col_character(),
  demo_rlgn = col_character(),
  demo_sex = col_character()
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_character(),
  X1 = col_double(),
  epi_subj = col_double(),
  epi_demo_age = col_double(),
  epi_demo_ses_num = col_double(),
  epi_demo_howr_num = col_double(),
  epi_demo_affr_num = col_double(),
  epi_demo_tung_num = col_double(),
  epi_demo_ytng_num = col_double(),
  epi_demo_affr_cat = col_logical(),
  epi_demo_tung_cat = col_logical(),
  epi_demo_ytng_cat = col_logical(),
  epi_version = col_double(),
  epi_charc = col_double()
)
See spec(...) for full column specifications.
2 parsing failures.
 row       col expected  actual                           file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  p7_ctry = col_character(),
  p7_abs_check = col_character(),
  p7_dse_check = col_character(),
  p7_se_check = col_character(),
  p7_unev_check = col_character(),
  p7_exsen_check = col_character(),
  p7_por_check = col_character(),
  p7_mm_check = col_character(),
  p7_dem_sex = col_character(),
  p7_dem_pocc = col_character(),
  p7_dem_major = col_character(),
  p7_dem_ethnicity = col_character(),
  p7_dem_rur.urb = col_character(),
  p7_dem_affrd.basics = col_character(),
  p7_dem_religion = col_character(),
  p7_dem_church = col_character(),
  p7_dem_holy.tung.gif = col_character(),
  p7_abs_child.exp_cat = col_logical(),
  p7_abs_poetic_cat = col_logical(),
  p7_abs_tv.real_cat = col_logical()
  # ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")

More on recoding (“QTS2”)

d1r_byquestion <- read_csv("../data/study1r_byquestion.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  researcher = col_character(),
  country = col_character(),
  quad = col_character(),
  religion = col_character(),
  coder = col_character(),
  site = col_character()
)
See spec(...) for full column specifications.

d1r <- read_csv("../data/study1r.csv") %>%
  rename_at(vars(contains("score")), funs(gsub("score", "recoded", .))) %>%
  full_join(d1) %>%
  mutate(country = factor(country, levels = levels_country),
         site = factor(site, levels = levels_site),
         religion = factor(religion, levels = levels_religion),
         researcher = factor(researcher, levels = levels_researcher))

Parsed with column specification:
cols(
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  coder = col_character(),
  subject_id = col_double(),
  spirit_score = col_double(),
  other_score = col_double(),
  study = col_character(),
  spirit_score_std = col_double(),
  other_score_std = col_double(),
  spirit_score_std2 = col_double(),
  other_score_std2 = col_double()
)
Joining, by = c("researcher", "country", "site", "religion", "subject_id", "study")
Column `researcher` joining character vector and factor, coercing into character vectorColumn `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector

contrasts(d1r$country) <- contrasts_country
contrasts(d1r$site) <- contrasts_site
contrasts(d1r$religion) <- contrasts_religion

Averages, by score

d1r %>% 
  select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
  gather(scale, score, -c(country, subject_id)) %>%
  filter(grepl("spirit", scale) | grepl("other", scale)) %>%
  mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                grepl("score", scale) ~ "original",
                                TRUE ~ NA_character_),
         score_type = factor(score_type, levels = c("original", "recoded")),
         scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
                           grepl("other", scale) ~ "other extraordinary experiences",
                           TRUE ~ NA_character_),
         scale = factor(scale, levels = c("spiritual experiences", 
                                          "other extraordinary experiences"))) %>%
  ggplot(aes(x = scale, y = score, group = score_type)) +
  # geom_point(aes(color = score_type), alpha = 0.2,
  #            position = position_jitterdodge(jitter.width = 0.1,
  #                                            dodge.width = 0.5,
  #                                            jitter.height = 0.04)) +
  geom_pointrange(data = . %>%
                    group_by(scale, score_type) %>%
                    langcog::multi_boot_standard(col = "score", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
                  position = position_dodge(width = 0.5),
                  shape = 21, fatten = 6) +
  labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
       x = "Scale", y = "Score (range: 0-1)",
       fill = "Score type")

lmer(score ~ score_type * scale + (1 | subject_id), 
   data = d1r %>% 
     select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
     gather(scale, score, -c(country, subject_id)) %>%
     filter(grepl("spirit", scale) | grepl("other", scale)) %>%
     mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                   grepl("score", scale) ~ "original",
                                   TRUE ~ NA_character_),
            score_type = factor(score_type, levels = c("original", "recoded")),
            scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
                              grepl("other", scale) ~ "other extraordinary experiences",
                              TRUE ~ NA_character_),
            scale = factor(scale, levels = c("spiritual experiences", 
                                             "other extraordinary experiences"))),
   contrasts = list(score_type = "contr.sum", scale = "contr.sum")) %>%
  summary()

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: score ~ score_type * scale + (1 | subject_id)
   Data: 
d1r %>% select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%  
    gather(scale, score, -c(country, subject_id)) %>% filter(grepl("spirit",  
    scale) | grepl("other", scale)) %>% mutate(score_type = case_when(grepl("recoded",  
    scale) ~ "recoded", grepl("score", scale) ~ "original", TRUE ~  
    NA_character_), score_type = factor(score_type, levels = c("original",  
    "recoded")), scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",  
    grepl("other", scale) ~ "other extraordinary experiences",  
    TRUE ~ NA_character_), scale = factor(scale, levels = c("spiritual experiences",  
    "other extraordinary experiences")))

REML criterion at convergence: 362.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.37869 -0.68474  0.01642  0.55450  2.82615 

Random effects:
 Groups     Name        Variance Std.Dev.
 subject_id (Intercept) 0.04320  0.2078  
 Residual               0.05204  0.2281  
Number of obs: 1336, groups:  subject_id, 341

Fixed effects:
                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)         4.251e-01  1.289e-02  3.416e+02  32.992  < 2e-16 ***
score_type1        -4.429e-03  6.269e-03  1.003e+03  -0.706 0.480112    
scale1              2.433e-02  6.252e-03  9.972e+02   3.892 0.000106 ***
score_type1:scale1  6.545e-03  6.246e-03  9.954e+02   1.048 0.294921    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scr_t1 scale1
score_type1  0.004              
scale1      -0.006  0.006       
scr_typ1:s1  0.003 -0.010  0.003

d1r %>% 
  select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
  gather(scale, score, -c(country, subject_id)) %>%
  filter(grepl("spirit", scale) | grepl("other", scale)) %>%
  mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                grepl("score", scale) ~ "original",
                                TRUE ~ NA_character_),
         score_type = factor(score_type, levels = c("original", "recoded")),
         scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
                           grepl("other", scale) ~ "other exp.",
                           TRUE ~ NA_character_),
         scale = factor(scale, levels = c("spiritual exp.", 
                                          "other exp."))) %>%
  ggplot(aes(x = scale, y = score, group = score_type)) +
  facet_grid(~ country) +
  # geom_point(aes(color = score_type), alpha = 0.2,
  #            position = position_jitterdodge(jitter.width = 0.1,
  #                                            dodge.width = 0.5,
  #                                            jitter.height = 0.04)) +
  geom_pointrange(data = . %>%
                    group_by(scale, score_type, country) %>%
                    langcog::multi_boot_standard(col = "score", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
                  position = position_dodge(width = 0.5),
                  shape = 21, fatten = 6) +
  labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
       x = "Scale", y = "Score (range: 0-1)",
       fill = "Score type") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

d1r %>% 
  select(country, site, religion, subject_id, 
         ends_with("recoded"), ends_with("score")) %>%
  gather(scale, score, -c(country, site, religion, subject_id)) %>%
  filter(grepl("spirit", scale) | grepl("other", scale)) %>%
  mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                grepl("score", scale) ~ "original",
                                TRUE ~ NA_character_),
         score_type = factor(score_type, levels = c("original", "recoded")),
         scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
                           grepl("other", scale) ~ "other exp.",
                           TRUE ~ NA_character_),
         scale = factor(scale, levels = c("spiritual exp.", 
                                          "other exp."))) %>%
  ggplot(aes(x = scale, y = score, group = score_type)) +
  facet_grid(rows = vars(site, religion), cols = vars(country)) +
  geom_point(aes(color = score_type), alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.5,
                                             jitter.height = 0.04)) +
  geom_pointrange(data = . %>%
                    group_by(scale, score_type, country, site, religion) %>%
                    langcog::multi_boot_standard(col = "score", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
                  position = position_dodge(width = 0.5),
                  shape = 21, fatten = 6) +
  labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
       x = "Scale", y = "Score (range: 0-1)",
       fill = "Score type", color = "Score type") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

d1r %>% 
  select(country, subject_id, coder, ends_with("recoded")) %>%
  gather(scale, score, -c(country, subject_id, coder)) %>%
  mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                grepl("score", scale) ~ "original",
                                TRUE ~ NA_character_),
         score_type = factor(score_type, levels = c("original", "recoded")),
         scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
                           grepl("other", scale) ~ "other exp.",
                           TRUE ~ NA_character_),
         scale = factor(scale, levels = c("spiritual exp.", 
                                          "other exp."))) %>%
  filter(grepl("spirit", scale) | grepl("other", scale)) %>%
  ggplot(aes(x = scale, y = score, group = coder)) +
  facet_grid(~ country) +
  geom_point(aes(color = coder), alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.5,
                                             jitter.height = 0.04)) +
  geom_pointrange(data = . %>%
                    group_by(scale, coder, country) %>%
                    langcog::multi_boot_standard(col = "score", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = coder),
                  shape = 21, fatten = 5,
                  position = position_dodge(width = 0.5)) +
  labs(title = "Comparing coders in recoded dataset (QTS2)",
       x = "Scale", y = "Score (range: 0-1)", fill = "Coder", color = "Coder") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

lmer(score ~ scale * coder + (1 | subject_id), 
   data = d1r %>% 
     select(country, coder, subject_id, ends_with("recoded")) %>%
     distinct() %>%
     gather(scale, score, -c(country, subject_id, coder)) %>%
     filter(grepl("spirit", scale) | grepl("other", scale)) %>%
     mutate(scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
                              grepl("other", scale) ~ "other extraordinary experiences",
                              TRUE ~ NA_character_),
            scale = factor(scale, levels = c("spiritual experiences", 
                                             "other extraordinary experiences")),
            coder = factor(coder, levels = c("Nikki", "Maria", "Lucy"))),
   contrasts = list(scale = "contr.sum", 
                    coder = cbind("_NvGM" = c(1, 0, -1),
                                  "_MvGM" = c(0, 1, -1)))) %>%
  summary()

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: score ~ scale * coder + (1 | subject_id)
   Data: d1r %>% select(country, coder, subject_id, ends_with("recoded")) %>%  
    distinct() %>% gather(scale, score, -c(country, subject_id,  
    coder)) %>% filter(grepl("spirit", scale) | grepl("other",  
    scale)) %>% mutate(scale = case_when(grepl("spirit", scale) ~  
    "spiritual experiences", grepl("other", scale) ~ "other extraordinary experiences",  
    TRUE ~ NA_character_), scale = factor(scale, levels = c("spiritual experiences",  
    "other extraordinary experiences")), coder = factor(coder,  
    levels = c("Nikki", "Maria", "Lucy")))

REML criterion at convergence: 328.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.76292 -0.68127  0.03902  0.47704  2.27078 

Random effects:
 Groups     Name        Variance Std.Dev.
 subject_id (Intercept) 0.02397  0.1548  
 Residual               0.06998  0.2645  
Number of obs: 671, groups:  subject_id, 340

Fixed effects:
                   Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)         0.42875    0.01515 341.70286  28.298  < 2e-16 ***
scale1              0.02745    0.01174 339.27621   2.338  0.01995 *  
coder_NvGM          0.02927    0.01839 340.24072   1.591  0.11246    
coder_MvGM         -0.06042    0.02038 340.51182  -2.965  0.00324 ** 
scale1:coder_NvGM  -0.03378    0.01423 337.85385  -2.374  0.01818 *  
scale1:coder_MvGM   0.01857    0.01577 338.11808   1.178  0.23978    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scale1 cd_NGM cd_MGM s1:_NG
scale1      -0.022                            
coder_NvGM  -0.434  0.015                     
coder_MvGM  -0.142  0.008 -0.173              
scl1:cd_NGM  0.015 -0.437 -0.017 -0.004       
scl1:cd_MGM  0.008 -0.145 -0.004 -0.018 -0.169

d1r %>% 
  select(country, site, religion, subject_id, coder, ends_with("recoded")) %>%
  gather(scale, score, -c(country, site, religion, subject_id, coder)) %>%
  mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
                                grepl("score", scale) ~ "original",
                                TRUE ~ NA_character_),
         score_type = factor(score_type, levels = c("original", "recoded")),
         scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
                           grepl("other", scale) ~ "other exp.",
                           TRUE ~ NA_character_),
         scale = factor(scale, levels = c("spiritual exp.", 
                                          "other exp."))) %>%
  filter(grepl("spirit", scale) | grepl("other", scale)) %>%
  ggplot(aes(x = scale, y = score, group = coder)) +
  facet_grid(rows = vars(site, religion), cols = vars(country)) +
  geom_point(aes(color = coder), alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1,
                                             dodge.width = 0.5,
                                             jitter.height = 0.04)) +
  geom_pointrange(data = . %>%
                    group_by(scale, coder, country, site, religion) %>%
                    langcog::multi_boot_standard(col = "score", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = coder),
                  shape = 21, fatten = 5,
                  position = position_dodge(width = 0.5)) +
  geom_text(data = . %>% count(country, site, religion, coder, scale),
             aes(y = 0, label = paste0("n=", n), color = coder), 
            position = position_dodge(width = 1)) +
  labs(title = "Comparing coders in recoded dataset (QTS2)",
       x = "Scale", y = "Score (range: 0-1)", fill = "Coder", color = "Coder") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

Averages, by question

d1r_byquestion %>% 
  select(country, site, religion, subject_id, 
         !!s1_var_spirit, !!s1_var_other) %>%
  gather(question, recoded, -c(country, site, religion, subject_id)) %>%
  full_join(d1_byquestion %>%
              select(country, site, religion, subject_id, 
                     !!s1_var_spirit, !!s1_var_other) %>%
              gather(question, original, -c(country, site, religion, subject_id))) %>%
  gather(response_type, response, c(original, recoded)) %>%
  mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other)),
         response_type = factor(response_type, levels = c("original", "recoded"))) %>%
  ggplot(aes(x = question, y = response, group = response_type)) +
  geom_pointrange(data = . %>%
                    group_by(question, response_type) %>%
                    langcog::multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper,
                      fill = response_type),
                  position = position_dodge(width = 0.5),
                  shape = 21, fatten = 6) +
  labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
       x = "Question", y = "Response (range: 0-1)",
       fill = "Response type") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector

d1r_byquestion %>% 
  select(country, site, religion, subject_id, 
         !!s1_var_spirit, !!s1_var_other) %>%
  gather(question, recoded, -c(country, site, religion, subject_id)) %>%
  full_join(d1_byquestion %>%
              select(country, site, religion, subject_id, 
                     !!s1_var_spirit, !!s1_var_other) %>%
              gather(question, original, -c(country, site, religion, subject_id))) %>%
  gather(response_type, response, c(original, recoded)) %>%
  mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other)),
         response_type = factor(response_type, levels = c("original", "recoded"))) %>%
  filter(question == "seehearnotgod") %>%
  count(response_type, response) %>%
  spread(response, n)

Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector

d1r_byquestion %>% 
  select(country, site, religion, subject_id, 
         !!s1_var_spirit, !!s1_var_other) %>%
  gather(question, recoded, -c(country, site, religion, subject_id)) %>%
  full_join(d1_byquestion %>%
              select(country, site, religion, subject_id, 
                     !!s1_var_spirit, !!s1_var_other) %>%
              gather(question, original, -c(country, site, religion, subject_id))) %>%
  mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other))) %>%
  filter(question == "seehearnotgod") %>%
  mutate(match = (original == recoded)) %>%
  count(match) %>%
  mutate(prop = n/sum(n))

Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector

Comparing first question to interviewer judgments when available

d1_check_questions <- c("godvoxaloud", "godviavisions", "godviatouch", 
                        "godviasmell", "neartangiblegod", "presencenotgod", 
                        "presencedemon", "beingentbody", "seehearnotgod", 
                        "whitelight", "spiritbeingencounter", "voxwhenalone", 
                        "seethingscornereye")

d1_check <- d1_byquestion %>% 
  select(country, subject_id, ends_with("_judge")) %>%
  gather(question, judgment, -c(country, subject_id)) %>%
  mutate(question = gsub("_judge", "", question)) %>%
  full_join(d1_byquestion %>% 
              select(country, subject_id, !!d1_check_questions) %>%
              gather(question, response, -c(country, subject_id)))

Joining, by = c("country", "subject_id", "question")

d1_check %>%
  filter(question %in% d1_check_questions) %>%
  mutate(agree = (judgment == response),
         agree = ifelse(is.na(agree), 0, agree)) %>%
  group_by(question) %>%
  summarise(percent_agree = mean(agree, na.rm = T)) %>%
  arrange(percent_agree)

s1_check_icc_df <- data.frame(question = character(), icc = numeric())
for(i in d1_check_questions){
  temp_icc <- icc_fun(d1_check, var_name = i, 
                      var1 = "judgment", var2 = "response")
  temp_res <- data.frame(question = i, icc = temp_icc)
  s1_check_icc_df <- full_join(s1_check_icc_df, temp_res)
  rm(temp_res, i)
}

Joining, by = c("question", "icc")
Column `question` joining factors with different levels, coercing to character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vector

s1_check_icc_df <- s1_check_icc_df %>%
  mutate(question = factor(question, levels = d1_check_questions),
         reliability = case_when(icc < 0.5 ~ "poor",
                                 icc < 0.75 ~ "moderate",
                                 icc < 0.9 ~ "good", 
                                 icc <= 1 ~ "excellent",
                                 TRUE ~ NA_character_),
         reliability = factor(reliability, 
                              levels = c("poor", "moderate", 
                                         "good", "excellent")),
         scale = case_when(question %in% s1_var_spirit ~ "spiritual experiences",
                           question %in% s1_var_other ~ "other experiences",
                           TRUE ~ NA_character_),
         scale = factor(scale, 
                        levels = c("spiritual experiences", "other experiences")))

s1_check_icc_df %>% arrange(icc)

s1_check_icc_df %>% 
  mutate(scale = recode_factor(scale,
                               "spiritual experiences" = "spiritual",
                               "other experiences" = "other")) %>%
  ggplot(aes(x = question, y = icc, color = reliability)) + 
  facet_grid(~ scale, scales = "free", space = "free") +
  geom_point() +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2)) +
  scale_color_brewer(palette = "RdYlGn", direction = 1, drop = F) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1,
                                   color = ifelse(s1_check_icc_df$reliability == "poor",
                                                  "#d7191c", "black"))) +
  labs(title = "Comparing initial responses vs. later 'judgments'",
       x = "Question", y = "Intraclass correlation coefficient (ICC)")

Percent with particular experiences

d1_byquestion %>% 
  count(godvoxaloud) %>% 
  mutate(godvoxaloud = recode_factor(godvoxaloud, 
                                     "0" = "no", "0.5" = "maybe", "1" = "yes"),
         proportion = round(n/sum(n), 2))

d1_byquestion %>% 
  filter(!is.na(godvoxaloud)) %>%
  count(country, godvoxaloud) %>% 
  mutate(godvoxaloud = recode_factor(godvoxaloud, 
                                     "0" = "no", "0.5" = "maybe", "1" = "yes")) %>%
  complete(country, nesting(godvoxaloud), fill = list(n = 0)) %>%
  spread(godvoxaloud, n) %>%
  mutate(total = no + maybe + yes) %>%
  janitor::adorn_totals() %>%
  mutate_at(vars(no, maybe, yes), funs(round(./total, 2)*100)) %>%
  select(-total) %>%
  data.frame() %>%
  mutate_at(vars(no, maybe, yes), funs(paste0(., "%"))) %>%
  kable() %>%
  kable_styling() %>%
  row_spec(6, bold = T)

Column `country` has different attributes on LHS and RHS of join

country	no	maybe	yes
US	86%	1%	13%
Ghana	33%	0%	67%
Thailand	58%	2%	40%
China	50%	3%	47%
Vanuatu	54%	9%	37%
Total	57%	3%	41%

d1_byquestion %>% 
  filter(!is.na(godviavisions)) %>%
  count(country, godviavisions) %>% 
  mutate(godviavisions = recode_factor(godviavisions, 
                                       "0" = "no", "0.5" = "maybe", "1" = "yes")) %>%
  complete(country, nesting(godviavisions), fill = list(n = 0)) %>%
  spread(godviavisions, n) %>%
  mutate(total = no + maybe + yes) %>%
  janitor::adorn_totals() %>%
  mutate_at(vars(no, maybe, yes), funs(round(./total, 2)*100)) %>%
  select(-total) %>%
  data.frame() %>%
  mutate_at(vars(no, maybe, yes), funs(paste0(., "%"))) %>%
  kable() %>%
  kable_styling() %>%
  row_spec(6, bold = T)

Column `country` has different attributes on LHS and RHS of join

country	no	maybe	yes
US	78%	0%	22%
Ghana	41%	0%	59%
Thailand	55%	2%	43%
China	78%	0%	22%
Vanuatu	61%	0%	39%
Total	63%	0%	37%

TML’s additional analyses

Focus on godvoxaloud and godviavisions

d1_tml1 <- d1 %>% select(country, subject_id, abs_score, por_score) %>%
  full_join(d1_byquestion %>% 
              select(country, subject_id, godvoxaloud, godviavisions) %>%
              gather(question, response, -c(country, subject_id)) %>%
              mutate(response_yn = case_when(response < 1 ~ 0,
                                             response == 1 ~ 1,
                                             TRUE ~ NA_real_)))

Joining, by = c("country", "subject_id")

glmer(response_yn ~ scale(por_score) + (1 | question) + (1 | country), 
      family = "binomial", d1_tml1) %>% summary()

boundary (singular) fit: see ?isSingular

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(por_score) + (1 | question) + (1 | country)
   Data: d1_tml1

     AIC      BIC   logLik deviance df.resid 
   769.4    787.2   -380.7    761.4      620 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7952 -0.7330 -0.4624  0.8706  2.4182 

Random effects:
 Groups   Name        Variance Std.Dev.
 country  (Intercept) 0.297    0.545   
 question (Intercept) 0.000    0.000   
Number of obs: 624, groups:  country, 5; question, 2

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)   
(Intercept)       -0.5236     0.2595  -2.018  0.04360 * 
scale(por_score)   0.3109     0.1090   2.851  0.00436 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(pr_scr) -0.006
convergence code: 0
boundary (singular) fit: see ?isSingular

glmer(response_yn ~ scale(abs_score) + (1 | question) + (1 | country), 
      family = "binomial", d1_tml1) %>% summary()

boundary (singular) fit: see ?isSingular

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | question) + (1 | country)
   Data: d1_tml1

     AIC      BIC   logLik deviance df.resid 
   745.3    762.9   -368.6    737.3      603 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6946 -0.7569 -0.4474  0.9293  2.7417 

Random effects:
 Groups   Name        Variance Std.Dev.
 country  (Intercept) 0.4476   0.669   
 question (Intercept) 0.0000   0.000   
Number of obs: 607, groups:  country, 5; question, 2

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.53208    0.31285  -1.701 0.088989 .  
scale(abs_score)  0.33016    0.09554   3.456 0.000549 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(bs_scr) -0.028
convergence code: 0
boundary (singular) fit: see ?isSingular

glmer(response_yn ~ scale(por_score) + (1 | country), 
      family = "binomial", d1_tml1 %>% filter(question == "godvoxaloud")) %>% 
  summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(por_score) + (1 | country)
   Data: d1_tml1 %>% filter(question == "godvoxaloud")

     AIC      BIC   logLik deviance df.resid 
   385.8    397.0   -189.9    379.8      310 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7575 -0.7684 -0.4053  0.7893  2.5614 

Random effects:
 Groups  Name        Variance Std.Dev.
 country (Intercept) 0.585    0.7649  
Number of obs: 313, groups:  country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)       -0.4694     0.3651  -1.286    0.199
scale(por_score)   0.2172     0.1551   1.401    0.161

Correlation of Fixed Effects:
            (Intr)
scl(pr_scr) 0.000 

glmer(response_yn ~ scale(abs_score) + (1 | country), 
      family = "binomial", d1_tml1 %>% filter(question == "godvoxaloud")) %>% 
  summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | country)
   Data: d1_tml1 %>% filter(question == "godvoxaloud")

     AIC      BIC   logLik deviance df.resid 
   375.5    386.7   -184.8    369.5      302 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6176 -0.7766 -0.3915  0.8313  2.8594 

Random effects:
 Groups  Name        Variance Std.Dev.
 country (Intercept) 0.7087   0.8418  
Number of obs: 305, groups:  country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)       -0.4886     0.3986  -1.226    0.220
scale(abs_score)   0.2052     0.1343   1.529    0.126

Correlation of Fixed Effects:
            (Intr)
scl(bs_scr) -0.019

glmer(response_yn ~ scale(por_score) + (1 | country), 
      family = "binomial", d1_tml1 %>% filter(question == "godviavisions")) %>% 
  summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(por_score) + (1 | country)
   Data: d1_tml1 %>% filter(question == "godviavisions")

     AIC      BIC   logLik deviance df.resid 
   384.3    395.5   -189.1    378.3      308 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6899 -0.7062 -0.5004  0.9736  2.2955 

Random effects:
 Groups  Name        Variance Std.Dev.
 country (Intercept) 0.09818  0.3133  
Number of obs: 311, groups:  country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)   
(Intercept)       -0.6004     0.1891  -3.175  0.00150 **
scale(por_score)   0.5138     0.1657   3.100  0.00194 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(pr_scr) -0.017

glmer(response_yn ~ scale(abs_score) + (1 | country), 
      family = "binomial", d1_tml1 %>% filter(question == "godviavisions")) %>% 
  summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | country)
   Data: d1_tml1 %>% filter(question == "godviavisions")

     AIC      BIC   logLik deviance df.resid 
   371.7    382.9   -182.9    365.7      299 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6618 -0.7152 -0.4837  0.9376  2.9033 

Random effects:
 Groups  Name        Variance Std.Dev.
 country (Intercept) 0.3431   0.5857  
Number of obs: 302, groups:  country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -0.6102     0.2933  -2.081 0.037477 *  
scale(abs_score)   0.4761     0.1387   3.432 0.000599 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(bs_scr) -0.066

Focus on other presence questions

d1_tml2 <- d1 %>% select(country, subject_id, abs_score, por_score) %>%
  full_join(d1_byquestion %>% 
              select(country, subject_id, neartangiblegod, presencenotgod, 
                     presencedemon, beingentbody, seehearnotgod, 
                     spiritbeingencounter) %>%
              gather(question, response, -c(country, subject_id)) %>%
              mutate(response_yn = case_when(response < 1 ~ 0,
                                             response == 1 ~ 1,
                                             TRUE ~ NA_real_)))

Joining, by = c("country", "subject_id")

glmer(response_yn ~ scale(por_score) + (1 | question) + (1 | country), 
      family = "binomial", d1_tml2) %>% summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(por_score) + (1 | question) + (1 | country)
   Data: d1_tml2

     AIC      BIC   logLik deviance df.resid 
  2025.9   2048.1  -1009.0   2017.9     1875 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7229 -0.6165 -0.3640  0.7589  4.0868 

Random effects:
 Groups   Name        Variance Std.Dev.
 question (Intercept) 0.4731   0.6878  
 country  (Intercept) 0.6302   0.7939  
Number of obs: 1879, groups:  question, 6; country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.75940    0.45630  -1.664   0.0961 .  
scale(por_score)  0.44752    0.06785   6.596 4.22e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(pr_scr) -0.012

glmer(response_yn ~ scale(abs_score) + (1 | question) + (1 | country), 
      family = "binomial", d1_tml2) %>% summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | question) + (1 | country)
   Data: d1_tml2

     AIC      BIC   logLik deviance df.resid 
  1978.2   2000.3   -985.1   1970.2     1822 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6347 -0.6291 -0.3769  0.7575  4.4684 

Random effects:
 Groups   Name        Variance Std.Dev.
 question (Intercept) 0.4829   0.6949  
 country  (Intercept) 0.6679   0.8172  
Number of obs: 1826, groups:  question, 6; country, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.79078    0.46635  -1.696   0.0899 .  
scale(abs_score)  0.28715    0.05964   4.814 1.48e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scl(bs_scr) -0.015

Cultural invitation vs. individual difference

Amount of variability each predictor, within each site

d1_variance <- d1 %>% 
  select(country, site, religion, subject_id, abs_score, por_score) %>%
  gather(var, score, ends_with("_score")) %>%
  filter(!is.na(score)) %>%
  mutate(score_rescaled = case_when(var == "abs_score" ~ score/1,
                                    var == "por_score" ~ score/3)) # called pv_score in study 4

bartlett.test(score_rescaled ~ var, 
              data = d1_variance %>% filter(country == "US"))

    Bartlett test of homogeneity of variances

data:  score_rescaled by var
Bartlett's K-squared = 2.3384, df = 1, p-value = 0.1262

d1_variance %>% 
  filter(country == "US") %>%
  group_by(var) %>%
  summarise(variance = var(score_rescaled))

bartlett.test(score_rescaled ~ var, 
              data = d1_variance %>% filter(country == "Ghana"))

    Bartlett test of homogeneity of variances

data:  score_rescaled by var
Bartlett's K-squared = 4.6186, df = 1, p-value = 0.03163

d1_variance %>% 
  filter(country == "Ghana") %>%
  group_by(var) %>%
  summarise(variance = var(score_rescaled))

bartlett.test(score_rescaled ~ var, 
              data = d1_variance %>% filter(country == "Thailand"))

    Bartlett test of homogeneity of variances

data:  score_rescaled by var
Bartlett's K-squared = 0.8459, df = 1, p-value = 0.3577

d1_variance %>% 
  filter(country == "Thailand") %>%
  group_by(var) %>%
  summarise(variance = var(score_rescaled))

bartlett.test(score_rescaled ~ var, 
              data = d1_variance %>% filter(country == "China"))

    Bartlett test of homogeneity of variances

data:  score_rescaled by var
Bartlett's K-squared = 9.3242, df = 1, p-value = 0.002262

d1_variance %>% 
  filter(country == "China") %>%
  group_by(var) %>%
  summarise(variance = var(score_rescaled))

bartlett.test(score_rescaled ~ var, 
              data = d1_variance %>% filter(country == "Vanuatu"))

    Bartlett test of homogeneity of variances

data:  score_rescaled by var
Bartlett's K-squared = 16.04, df = 1, p-value = 6.202e-05

d1_variance %>% 
  filter(country == "Vanuatu") %>%
  group_by(var) %>%
  summarise(variance = var(score_rescaled))

d1_variance %>% 
  group_by(var, country) %>%
  summarise(variance = var(score_rescaled)) %>%
  ungroup() %>%
  spread(var, variance) %>%
  rename(Country = country, Absorption = abs_score, 
         `Porosity Vignettes` = por_score) %>%
  kable(digits = 2) %>%
  kable_styling()

Country	Absorption	Porosity Vignettes
US	0.04	0.03
Ghana	0.05	0.03
Thailand	0.03	0.02
China	0.05	0.02
Vanuatu	0.03	0.01

d1_variance %>%
  group_by(var, country) %>%
  mutate(score_cent = scale(score_rescaled, scale = F)) %>%
  ungroup() %>%
  mutate(var = recode_factor(var,
                             "por_score" = "Porosity Vignettes",
                             "abs_score" = "Absorption")) %>%
  ggplot(aes(x = score_cent, fill = var, color = var, lty = var)) +
  facet_wrap(country ~ ., ncol = 3) +
  geom_density(alpha = 0.1) +
  scale_linetype_manual(values = c(1, 2)) +
  scale_color_brewer(palette = "Set1", direction = -1) +
  scale_fill_brewer(palette = "Set1", direction = -1) +
  labs(x = "Score (rescaled to run from 0-1, centered at the mean by country)",
       color = "Scale", fill = "Scale", lty = "Scale") +
  theme(legend.position = "bottom")

d1_variance %>%
  group_by(var, country, site, religion) %>%
  mutate(score_cent = scale(score_rescaled, scale = F)) %>%
  ungroup() %>%
  mutate(var = recode_factor(var,
                             "por_score" = "Porosity Vignettes",
                             "abs_score" = "Absorption")) %>%
  ggplot(aes(x = score_cent, fill = var, color = var, lty = var)) +
  facet_grid(cols = vars(country), rows = vars(site, religion)) +
  geom_density(alpha = 0.1) +
  scale_linetype_manual(values = c(1, 2)) +
  scale_color_brewer(palette = "Set1", direction = -1) +
  scale_fill_brewer(palette = "Set1", direction = -1) +
  labs(x = "Score (rescaled to run from 0-1, centered at the mean by country)",
       color = "Scale", fill = "Scale", lty = "Scale") +
  theme(legend.position = "bottom")

Variance accounted for by country

r7 <- lm(por_score_std ~ country, d1)
regtab_fun(r7, country_var1 = "country_gh", country_var2 = "country_th", country_var3 = "country_ch", country_var4 = "country_vt")

rsquared(r7)

r9 <- lm(abs_score_std ~ country, d1)
rsquared(r9)